<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"
                "http://www.w3.org/TR/REC-html40/loose.dtd">
<html>
<head>
  <title>Description of CellsortPlotPCspectrum</title>
  <meta name="keywords" content="CellsortPlotPCspectrum">
  <meta name="description" content="function CellsortPlotPCspectrum(fn, CovEvals, pcuse)">
  <meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
  <meta name="generator" content="m2html &copy; 2003 Guillaume Flandin">
  <meta name="robots" content="index, follow">
  <link type="text/css" rel="stylesheet" href="m2html.css">
</head>
<body>
<a name="_top"></a>
<div> <a href="index.html">CellSort 1.0</a> &gt; CellsortPlotPCspectrum.m</div>

<!--<table width="100%"><tr><td align="left"><a href="../index.html"><img alt="<" border="0" src="images/left.png">&nbsp;Master index</a></td>
<td align="right"><a href="index.html">Index for CellSort 1.0&nbsp;<img alt=">" border="0" src="images/right.png"></a></td></tr></table>-->

<h1>CellsortPlotPCspectrum
</h1>

<h2><a name="_name"></a>PURPOSE <a href="#_top"><img alt="^" border="0" src="images/up.png"></a></h2>
<div class="box"><strong>function CellsortPlotPCspectrum(fn, CovEvals, pcuse)</strong></div>

<h2><a name="_synopsis"></a>SYNOPSIS <a href="#_top"><img alt="^" border="0" src="images/up.png"></a></h2>
<div class="box"><strong>function CellsortPlotPCspectrum(fn, CovEvals, pcuse) </strong></div>

<h2><a name="_description"></a>DESCRIPTION <a href="#_top"><img alt="^" border="0" src="images/up.png"></a></h2>
<div class="fragment"><pre class="comment">CellsortPlotPCspectrum(fn, CovEvals, pcuse)

 Plot the principal component (PC) spectrum and compare with the
 corresponding random-matrix noise floor

 Inputs:
   fn - movie file name. Must be in TIFF format.
   CovEvals - eigenvalues of the covariance matrix
   pcuse - [optional] - indices of PCs included in dimensionally reduced
   data set

 Eran Mukamel, Axel Nimmerjahn and Mark Schnitzer, 2009
 Email: eran@post.harvard.edu, mschnitz@stanford.edu</pre></div>

<!-- crossreference -->
<h2><a name="_cross"></a>CROSS-REFERENCE INFORMATION <a href="#_top"><img alt="^" border="0" src="images/up.png"></a></h2>
This function calls:
<ul style="list-style-image:url(../matlabicon.gif)">
</ul>
This function is called by:
<ul style="list-style-image:url(../matlabicon.gif)">
</ul>
<!-- crossreference -->

<h2><a name="_subfunctions"></a>SUBFUNCTIONS <a href="#_top"><img alt="^" border="0" src="images/up.png"></a></h2>
<ul style="list-style-image:url(../matlabicon.gif)">
<li><a href="#_sub1" class="code">function formataxes</a></li><li><a href="#_sub2" class="code">function j = tiff_frames(fn)</a></li></ul>
<h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" src="images/up.png"></a></h2>
<div class="fragment"><pre>0001 <a name="_sub0" href="#_subfunctions" class="code">function CellsortPlotPCspectrum(fn, CovEvals, pcuse)</a>
0002 <span class="comment">% function CellsortPlotPCspectrum(fn, CovEvals, pcuse)</span>
0003 <span class="comment">%</span>
0004 <span class="comment">% Plot the principal component (PC) spectrum and compare with the</span>
0005 <span class="comment">% corresponding random-matrix noise floor</span>
0006 <span class="comment">%</span>
0007 <span class="comment">% Inputs:</span>
0008 <span class="comment">%   fn - movie file name. Must be in TIFF format.</span>
0009 <span class="comment">%   CovEvals - eigenvalues of the covariance matrix</span>
0010 <span class="comment">%   pcuse - [optional] - indices of PCs included in dimensionally reduced</span>
0011 <span class="comment">%   data set</span>
0012 <span class="comment">%</span>
0013 <span class="comment">% Eran Mukamel, Axel Nimmerjahn and Mark Schnitzer, 2009</span>
0014 <span class="comment">% Email: eran@post.harvard.edu, mschnitz@stanford.edu</span>
0015 
0016 <span class="keyword">if</span> nargin&lt;3
0017     pcuse = [];
0018 <span class="keyword">end</span>
0019 
0020 [pixw,pixh] = size(imread(fn,1));
0021 npix = pixw*pixh;
0022 nt = <a href="#_sub2" class="code" title="subfunction j = tiff_frames(fn)">tiff_frames</a>(fn);
0023 
0024 <span class="comment">% Random matrix prediction (Sengupta &amp; Mitra)</span>
0025 p1 = npix; <span class="comment">% Number of pixels</span>
0026 q1 = nt; <span class="comment">% Number of time frames</span>
0027 q = max(p1,q1);
0028 p = min(p1,q1);
0029 sigma = 1;
0030 lmax = sigma*sqrt(p+q + 2*sqrt(p*q));
0031 lmin = sigma*sqrt(p+q - 2*sqrt(p*q));
0032 lambda = [lmin: (lmax-lmin)/100.0123423421: lmax];
0033 rho = (1./(pi*lambda*(sigma^2))).*sqrt((lmax^2-lambda.^2).*(lambda.^2-lmin^2));
0034 rho(isnan(rho)) = 0;
0035 rhocdf = cumsum(rho)/sum(rho);
0036 noiseigs = interp1(rhocdf, lambda, [p:-1:1]'/p, <span class="string">'linear'</span>, <span class="string">'extrap'</span>).^2 ;
0037 
0038 <span class="comment">% Normalize the PC spectrum</span>
0039 normrank = min(nt-1,length(CovEvals));
0040 pca_norm = CovEvals*noiseigs(normrank) / (CovEvals(normrank)*noiseigs(1));
0041 
0042 clf
0043 plot(pca_norm, <span class="string">'o-'</span>, <span class="string">'Color'</span>, [1,1,1]*0.3, <span class="string">'MarkerFaceColor'</span>, [1,1,1]*0.3, <span class="string">'LineWidth'</span>,2)
0044 hold on
0045 plot(noiseigs / noiseigs(1), <span class="string">'b-'</span>, <span class="string">'LineWidth'</span>,2)
0046 plot(2*noiseigs / noiseigs(1), <span class="string">'b--'</span>, <span class="string">'LineWidth'</span>,2)
0047 <span class="keyword">if</span> ~isempty(pcuse)
0048     plot(pcuse, pca_norm(pcuse), <span class="string">'rs'</span>, <span class="string">'LineWidth'</span>,2)
0049 <span class="keyword">end</span>
0050 hold off
0051 <a href="#_sub1" class="code" title="subfunction formataxes">formataxes</a>
0052 set(gca,<span class="string">'XScale'</span>,<span class="string">'log'</span>,<span class="string">'YScale'</span>,<span class="string">'log'</span>, <span class="string">'Color'</span>,<span class="string">'none'</span>)
0053 xlabel(<span class="string">'PC rank'</span>)
0054 ylabel(<span class="string">'Normalized variance'</span>)
0055 axis tight
0056 <span class="keyword">if</span> isempty(pcuse)
0057     legend(<span class="string">'Data variance'</span>,<span class="string">'Noise floor'</span>,<span class="string">'2 x Noise floor'</span>)
0058 <span class="keyword">else</span>
0059     legend(<span class="string">'Data variance'</span>,<span class="string">'Noise floor'</span>,<span class="string">'2 x Noise floor'</span>,<span class="string">'Retained PCs'</span>)
0060 <span class="keyword">end</span>
0061 
0062 fntitle = fn;
0063 fntitle(fn==<span class="string">'_'</span>) = <span class="string">' '</span>;
0064 title(fntitle)
0065 
0066 <a name="_sub1" href="#_subfunctions" class="code">function formataxes</a>
0067 
0068 set(gca,<span class="string">'FontSize'</span>,12,<span class="string">'FontWeight'</span>,<span class="string">'bold'</span>,<span class="string">'FontName'</span>,<span class="string">'Helvetica'</span>,<span class="string">'LineWidth'</span>,2,<span class="string">'TickLength'</span>,[1,1]*.02,<span class="string">'tickdir'</span>,<span class="string">'out'</span>)
0069 set(gcf,<span class="string">'Color'</span>,<span class="string">'w'</span>,<span class="string">'PaperPositionMode'</span>,<span class="string">'auto'</span>)
0070 
0071 <a name="_sub2" href="#_subfunctions" class="code">function j = tiff_frames(fn)</a>
0072 <span class="comment">%</span>
0073 <span class="comment">% n = tiff_frames(filename)</span>
0074 <span class="comment">%</span>
0075 <span class="comment">% Returns the number of slices in a TIFF stack.</span>
0076 <span class="comment">%</span>
0077 <span class="comment">%</span>
0078 
0079 status = 1; j=0;
0080 jstep = 10^3;
0081 <span class="keyword">while</span> status
0082     <span class="keyword">try</span>
0083         j=j+jstep;
0084         imread(fn,j);
0085     <span class="keyword">catch</span>
0086         <span class="keyword">if</span> jstep&gt;1
0087             j=j-jstep;
0088             jstep = jstep/10;
0089         <span class="keyword">else</span>
0090             j=j-1;
0091             status = 0;
0092         <span class="keyword">end</span>
0093     <span class="keyword">end</span>
0094 <span class="keyword">end</span></pre></div>
<hr><address>Generated on Wed 29-Jul-2009 12:46:53 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/">m2html</a></strong> &copy; 2003</address>
</body>
</html>